unit TiUncertain;

interface

    Uses Math, Titrdefs, Solver;

{GasDev returns a random number from a normal dist. of given mean and StdDev}
Function GasDev(mean, dev : double) : double;

{Calcstats returns the mean, StdDev, skewness and Kurtosis in StatsOut}
Procedure CalcStats(Var X: MCConc; N : integer; Var StatsOut : Stats);

Procedure MonteMove( Var S : DataSet);       { Prepare arrays for MC }

Procedure OneConMonte(Var S : DataSet; Constraint : integer);

Procedure StoreMCFile(Var MCFullPoint : MCResult;       {Saves MC results}
                      FileName : string;                {for a single point}
                      ColMax, NumTrials : integer);

Procedure StoreMCSummary(Var SumStats : ConcStats;      {Saves Summary Stats}
                         FileName : string;             {For a single MC point}
                         NumConc: integer);

Function  MonteCarlo(Var S : DataSet; Var ThisPoint: MCResult; MaxReps : integer)         {Maximum # realizations}
                    : integer;               { Perform a single MC expt.}

Function CalcDistribution(Var ConcMat : MCResult;  Var Frequencies : MCFreqs;
       Var BinBounds : MCbounds; Conc : integer; NReps : integer;
                          NConc : integer) : Boolean;


Procedure SaveSensMat(Var SensMat : SensiMat; Filename : string; NConc : integer);       {number of Concentrations}

Procedure CalcSensMat(Var S: DataSet; Var SensMat: Sensimat; Delta : double);

Function  Predict(Order : integer) : Boolean; { Error prediction }
Function  MoSweep : integer;                 { Sweep with added MC }
Function  Surface : integer;                 { MC Titration for distribution}

{*****************************************************************************}

implementation

{ The procedures in the next block of code are all used to perform the }
{ Monte Carlo Simulations and sensitivity (derivative) calculations. }

{GasDev is a Gaussian random number generator from Numerical Recipes,
 page 203 (The Box-Muller transformation method).
 Returns a random number with mu = Mean and sigma = dev.}

 Function GasDev(mean, dev : double) : double;

 Var
   x, y, R, Factor : double;

 begin
   repeat
     x := 2*Random -1;
     y := 2*Random -1;
     R := x*x + y*y;
   until((R < 1) and (R <> 0));
   Factor := Sqrt(-2 * ln(R) / R);
   GasDev := Mean + Dev * x * Factor;
 end;                                      {End GasDev}

{GetConcs moves values of log[] from the component and species computational }
{arrays into the Target array}

Procedure GetConcs(Var Target: Concs);

Var com, spe : integer;

Begin
   For com := 1 to S.NComp do
      Target[com] := logX[com] / LogConvr;

   For spe := 1 to S.NSpec do begin
      If (conc[spe] > 0) then
        Target[spe + S.NComp] := ln(conc[spe])/ LogConvr
           else Target[spe + S.NComp] :=  -999;
     end;
end;

{ CalcStats calculates the mean, standard deviation, skewness and kurtosis of }
{ an array X with n elements. These values are returned in StatsOut }

Procedure CalcStats(Var X: MCConc; N : integer; Var StatsOut : Stats);

var    i            : integer;               {counter for array X}
       SumX, SumX2  : double;                {sums of X and X^2}
       Sum3, Sum4   : double;                {sums of (x-xbar)^3 and (x-xbar)^4}
       Xbar, SD     : double;                {Xbar = mean X}
       DeltaCheck   : double;                {check on SD calc}

Begin
  SumX := 0;                   {First set all sums to zero}
  SumX2 := 0;
  Sum3 := 0;
  Sum4 := 0;

  for i := 1 to n do
    begin
      SumX := SumX + X[i];             {Get sum of all x }
      SumX2 := SumX2 + X[i] * X[i]     {and all x^2 }
    end;

  DeltaCheck := SumX2 - (SumX * SumX) / n;
  If (DeltaCheck < 1e-14) then begin     {Deltacheck is very small, roundoff}
      StatsOut[Mean] := SumX / n;        {error may occur.  Avoid this by }
      StatsOut[StdDev] := 0;             {all higher moments to zero}
      StatsOut[Skewness] := 0;
      StatsOut[Kurtosis] := 0;
      Exit;
     end;

  SD := Sqrt(DeltaCheck / (n-1));        { calculate std. deviation}
  Xbar := SumX / n;                      { calculate mean }

   for i := 1 to N do                     {Now get 3rd and 4th moments of X}
     begin
       Sum3 := Sum3 + Power((X[i] - Xbar), 3);
       Sum4 := Sum4 + Power((X[i] - Xbar), 4);
     end;

   StatsOut[Mean] := Xbar;                {output statistics}
   StatsOut[StdDev] := SD;
   StatsOut[Skewness] := Sum3 / ((N-1) * SD * SD *SD);

   {Note: Kurtosis is 'normalized' to zero by subtracting 3 }
   StatsOut[Kurtosis] := (Sum4 / ((N-1) * SD * SD * SD * SD)) - 3;

End;                                      {end of CalcStats}

{ MonteMove places random values of log K, TotX, LogX into computational arrays}
{ for MC calculations}

Procedure MonteMove(Var S : DataSet);
Var
   com, spe : integer;           {Component and species calculators}

begin

  with S do begin                 { S is the data set }

    for com := 1 to NComp do begin

    { For all Known-total components, generate random Tot[] }
      if (Comps[com].CType = CTotDiss) then
        if (Comps[com].StdDev > 0) then
          Tot[com] := GasDev(Comps[com].TotalX, Comps[com].StdDev)
        else Tot[com] := Comps[com].TotalX;

     { For all known-activity components, generate random LogX[] }
      if (Comps[com].CType = CFixDiss) then
        if (Comps[com].StdDev > 0) then
         LogX[com] := GasDev(Comps[com].GuessX, Comps[com].StdDev) * LogConvr
        else LogX[com] := Comps[com].GuessX * LogConvr;

      end;                                       {end for ...}

    { For all species, generate random Log K }
    for spe := 1 to NSpec do begin
      If (Specs[spe].StdDev > 0) then
         LogKZero[spe] := GasDev(Specs[spe].LogK, Specs[spe].StdDev) * LogConvr
        else LogKZero[spe] := Specs[spe].LogK * LogConvr;
      end;

   end;     { End of the "With S " statement }
end;        { End of MonteMove }

{Procedure OneConMonte generates a random value for only one constraint, using}
{'mean' or 'true' values for all others.}

Procedure OneConMonte(Var S : DataSet; Constraint : integer);

Var
   spe : integer;                 {Species index}

begin
  InitArrays(False);              { initialize transfer true values to the }
                                  { computational arrays, retaining old freeX }

  with S do begin                              { S is the data set }
    If (Constraint > NComp) then begin         {Constraint is a log K}
        spe := Constraint - NComp;
        If (Specs[spe].StdDev > 0) then
         LogKZero[spe] := GasDev(Specs[spe].LogK, Specs[spe].StdDev) * LogConvr
       end                                     {end log K randomization}

     else begin                                {Constraint is a component}
       if ((Comps[Constraint].CType = CTotDiss) and     {Known Total component}
           (Comps[Constraint].StdDev > 0)) then
            Tot[Constraint] :=
            GasDev(Comps[Constraint].TotalX, Comps[Constraint].StdDev);

       if ((Comps[Constraint].CType = CFixDiss) and   {Free activity component}
           (Comps[Constraint].StdDev > 0)) then
            LogX[Constraint] :=
            GasDev(Comps[Constraint].GuessX, Comps[Constraint].StdDev) * LogConvr;

       end;                                     {end else}

   end;     { End of the "With S " statement }
end;        { End of OneConMonte }


{ StoreMCFile saves the Monte Carlo results matrix MCFullPoint (dimensions
{ [NumTrials] x [NumConc] as a text file called FileName. }

Procedure StoreMCFile(Var MCFullPoint : MCResult;
                      FileName : string;
                      ColMax, NumTrials : integer);

Var com, spe, mcn : integer;      {indices for components, species and trials}
    F             : TextFile;     {location to write to}

begin

    AssignFile(F, FileName);          { Open file for writing}
         {$I-}

    ReWrite(F);                       { start at the beginning of the file}

    Write(F, 'Trial ');

  {Write the names of the concentration variables calculated }

    For com := 1 to S.NComp do
        Write(F, S.Comps[com].Name,' ');
    For spe := 1 to S.NSpec do
        Write(F, S.Specs[spe].name, ' ');
    Writeln(F);

    {Now write the actual data, one line for each x value}

    For mcn := 1 to NumTrials do begin
       {str(, dummy);
       dummy := trimleft(dummy);  }
       Write(F, mcn, ' ');
       For spe := 1 to ColMax do
         Write(F, MCFullPoint[mcn, spe]:8:3, ' ');
       WriteLn(F);
      end;

     CloseFile(F);
   end;                           {end of StoreMCFile}

{StoreMCSummary saves the Monte Carlo summary stats for a single point }
{matrix SumStats, [NumConcs] x [4] in the text file filename)}

Procedure StoreMCSummary(Var SumStats : ConcStats;
                         FileName : string;
                         NumConc: integer);

Var com, spe   : integer;      {indices for components, species and trials}
    Stat       : StatType;     {index for statistics}
    F          : TextFile;     {location to write to}

begin

    AssignFile(F, FileName);          { Open file for writing}
         {$I-}

    ReWrite(F);                       { start at the beginning of the file}

    Writeln(F, 'Species TrueValue Mean StdDev Skewness Kurtosis');

  {Write the names of the concentration variables calculated }

    For com := 1 to S.NComp do begin
        Write(F, S.Comps[com].Name);
        For Stat := TrueValue to Kurtosis do
          Write(F, ' ', SumStats[com][Stat]:8:3);
        Writeln(F);
       end;                                    {end for com}
    For spe := 1 to S.NSpec do begin
        Write(F, S.Specs[spe].Name);
        For Stat := TrueValue to Kurtosis do
          Write(F, ' ', SumStats[Spe + S.NComp][Stat]:8:3);
        Writeln(F);
       end;

    CloseFile(F);
   end;                           {end of StoreMCFile}



{ MonteCarlo performs an error simulation by repeating the calculation  }
{ process with randomly generated values of TotX, FreeX and LogK, as    }
{ given by the mean and standard deviations in the data set S.          }
{ It calls InitArrays, MonteMove, Solidsolve                            }
{ Returns # failures to converge                                        }

Function MonteCarlo(Var S : DataSet;          {Input data}
                    Var ThisPoint: MCResult;  {Results matrix}
                    MaxReps : integer         {Maximum # realizations}
                    ): integer;               {returns 0 if successful}

Var
  mcn   : integer;               { current simulation # }
  com, spe  : integer;           { indices }
  Flag  : Boolean;               { Error flag for calculation }
  ErrorCount : integer;          { counts failures to converge }


begin
  InitArrays(TRUE);            { Initialize computational arrays}
  Flag := SolidSolve;          { Solve the system w/ 'true' constraints}

  With S do begin              { Save these as 'true' results }
      for com := 1 to NComp do          {For each Component, record log X}
         MCPointSum[com][TrueValue] := LogX[com]/logConvr;
      for spe := 1 to NSpec do begin    { For each species in the system S}
         If (conc[spe] > 0) then       { Record log C if it is present}
            MCPointSum[NComp + spe][TrueValue] := ln(conc[spe])/ LogConvr
                else MCPointSum[NComp +spe][TrueValue] := -999;
            end;                            { End for spe...}

       end;                    { end of with S}

  mcn := 0;                    { No realizations yet }
  ErrorCount := 0;             { No failures to converge }

  {Main loop for the MC simulation}
  Repeat
     mcn := mcn + 1;             { simulation counter }

     { Generate new random values }
     If OneCon then OneConMonte(S, ChosenCon)         {for one constraint}
        else MonteMove(S);                            {or all constraints}

     Flag := SolidSolve;         { Solve the present system }

     If (Flag = False) then begin          {Failed to converge, so}
        ErrorCount := ErrorCount + 1;      {increment the error counter}
        mcn := mcn-1;
      end

      else begin   {Read results from computational arrays into ThisPoint}
        With S do begin
          for com := 1 to NComp do          {For each Component, record log X}
              ThisPoint[mcn][com] := LogX[com]/logConvr;
          for spe := 1 to NSpec do begin    { For each species in the system S}
              If (conc[spe] > 0) then       { Record log C if it is present}
                 ThisPoint[mcn][NComp + spe] := ln(conc[spe])/ LogConvr
                    else Specs[spe].LogC := -999;
            end;                            { End for spe...}

       end;                                 {end with S}
      end;                                  {end of else...}

  Until((ErrorCount = ErrorMax) or (mcn = MaxReps));
  {End of repeat loop}

  MonteCarlo := ErrorCount;                 { Return # failures }
  MCDone := TRUE;

  end;                                      { End of MonteCarlo}

{Procedure SaveSensMat writes the sensitivity matrix SensMat into the text }
{file filename}

Procedure SaveSensMat(Var SensMat : SensiMat; {NConc x NConc Sensitivity matrix}
                      Filename : string;      {file to store in}
                      NConc : integer);       {number of Concentrations}

Var F : Text;
    com, spe, conc : integer;

begin
    AssignFile(F, FileName);          { Open file for writing}
         {$I-}

    ReWrite(F);                       { start at the beginning of the file}

    { Now write component and species names in the first row }
    Write(F, 'Conc. ');
    For com := 1 to S.NComp do
      Write(F, S.Comps[com].Name, ' ');
    For spe := 1 to S.NSpec do
      Write(F, S.Specs[spe].Name, ' ');
    WriteLn(F);

   {Write the values  of the first derivatives for all components}

    For com := 1 to S.NComp do begin
        Write(F, S.Comps[com].Name);
        For Conc := 1 to NConc do
          Write(F, ' ', SensMat[com][Conc]:8:3);
        Writeln(F);
       end;                                    {end for com loop}

    {Write the values  of the first derivatives for all components}
    For spe := 1 to S.NSpec do begin
        Write(F, S.Specs[spe].Name);
        For Conc := 1 to NConc do
          Write(F, ' ', SensMat[spe+S.NComp][Conc]:8:3);
        Writeln(F);
       end;

    CloseFile(F);

end;    {end of procedure StoreSensMat}

{CalcSensMat calculates a sensitivity matrix composed of the numerical first }
{derivatives of each concentration with respect to each constraint, }
{using values for each constraint X of X-delta and X + delta}
{where delta is the offset in log terms.}

Procedure CalcSensMat(Var S: DataSet;         {dataset to use}
                      Var SensMat: Sensimat;  {sensitivity matrix}
                      Delta : double);        {difference in constraint for }
                                              {calculating numerical derivative}

Var     com, spe, conc : integer; {counters/indices}
        NConc : integer;          {Number of concentrations, NConc + NSpec}
        DelFrac : double;         {difference in total X as a fraction }
        LowEnd, HighEnd : Concs;  {concentrations values for x-delta}
        Flag : Boolean;           {set true if delta X is zero}

begin
    {Initialize by calculating NConc and solving system once}
    NConc := S.NComp + S.NSpec;
    DelFrac := Power(10, delta);
    InitArrays(False);
    SolidSolve;

    {Now for each component, calculate X-delta, solve the system and save all }
    {concentrations.  Then calculate Derivative D = (C(X+delta)-C(X-delta)) /2}
    {and store in the sensitivity matrix}

    For com := 1 to S.NComp do begin
      {First get all concentrations for X-delta into Lowend}
      Initarrays(True);
      If (S.comps[com].CType = CTotDiss) then
          Tot[com] := S.comps[com].TotalX / DelFrac
         else LogX[com] := (S.Comps[com].GuessX-delta) * 2.3025851;
      SolidSolve;
      GetConcs(LowEnd);

      {Then get all concentrations for X+delta into highend}
      If (S.comps[com].CType = CTotDiss) then
         Tot[com] := S.comps[com].TotalX * DelFrac
        else LogX[com] := (S.Comps[com].GuessX+delta) * 2.3025851;
      SolidSolve;
      GetConcs(HighEnd);

      For conc := 1 to NConc do
         if ((S.comps[com].CType = CTotDiss) and (S.Comps[com].totalX = 0)) then
           SensMat[com][conc] := 999
           else SensMat[com][conc] := 0.5*(HighEnd[conc]-LowEnd[conc])/delta;

     end;

    For spe := 1 to S.NSpec do begin

      {First get all concentrations for X-delta into Lowend}
      Initarrays(True);
      LogKZero[spe] := (S.Specs[spe].LogK - delta) * 2.3025851;
      SolidSolve;
      GetConcs(LowEnd);

      {Then get all concentrations for X+delta into highend}
      LogKZero[spe] := (S.Specs[spe].LogK + delta) * 2.3025851;
      SolidSolve;
      GetConcs(HighEnd);

      For conc := 1 to NConc do
         SensMat[S.NComp + spe][conc] := 0.5*(HighEnd[conc]-LowEnd[conc])/delta;

     end;           {end of for spe loop}

end;

{ CalcDistribution determines the frequency distribution of the concentrations }
{ of species # conc in Monte Carlo Results matrix ConcMat (dimensions NReps by }
{ NConc), returning the values in the integer array frequencies.}

Function CalcDistribution(Var ConcMat : MCResult;       {MC single point result}
                          Var Frequencies : MCFreqs;    {counts/frequencies}
                          Var BinBounds : MCbounds;     {X values for bins}
                          Conc : integer;               {species of interest}
                          NReps : integer;              {number of MC trials}
                          NConc : integer) : Boolean;   {number of species}

Var   rep : integer;                {counts the number of trials}
      HighVal, LowVal : double;     {high and low concetration values }
      Range : double;               {range is high-low value}
      bin : integer;                {counts bins in histogram}
      LowBound : double;            {lowest value rounded down}

begin
  Result := TRUE;              {If nothing goes wrong, return true result}

  HighVal := -1000;            {set high bound very low}
  LowVal := 1000;              {and low bound very high}

  {For all trials, compare log [] with the bounds in order to find max & min}
  for rep := 1 to NReps do begin
      If (ConcMat[rep, Conc] > HighVal) then HighVal := ConcMat[rep, Conc];
      If (ConcMat[rep, Conc] < LowVal) then LowVal := ConcMat[rep, Conc];
    end;

  Range := HighVal - LowVal;

  {If all values are the same, no distribution is possible: Exit, return False}
  if (Range < 0.001) then begin
      Result := False;
      Exit;
     end;

  {Now calculate the bin boundaries for the histogram}

  BinWidth := (ceil(10*Range/MaxBins))/10;

  LowBound := floor(10*LowVal)/10;
  If ((Range/Binwidth) < 10) then BinWidth := Binwidth/2;
  If ((Range/Binwidth) < 20) then BinWidth := Binwidth/2;

  {Set Bin interval values in BinBounds and set all frequency counts to zero}
  For bin := 1 to MaxBins do begin
     BinBounds[bin] := LowBound + BinWidth * bin;
     Frequencies[bin] := 0;
    end;

  {Finally, count occurences to assign frequencies to each bin}
  For rep := 1 to NReps do begin
    bin := Floor((ConcMat[rep][conc] - LowBound)/ Binwidth);
    Frequencies[bin] := Frequencies[bin] + 1;
  end;

end;                                 {End of CalcDistribution}

{************************************************************************}
{                                                                        }
{   Error Prediction Code- predicts the standard deviations of           }
{                          Monte Carlo simulations                       }
{                                                                        }
{************************************************************************}

{ Function PREDICT-                                                      }
{ For a given data set, estimates the standard deviation of calculated   }
{ concentrations by calculating the numerical derivative dC/dX, dC/dT or }
{ dC/dK for each constraining variable (X, T or K). The squares of the   }
{ derivatives times the SD's of the constraints are summed in DeltaX and }
{ DeltaC, then the square root of the sum is calculated                  }

Function Predict(Order : integer) : Boolean;

                             { DMult is the fraction of the SD used in the }
                             { derivative routine. Order indicates the type }
                             { of derivative used }

Var OffX : Array[1..6] of Vector; { Log X values for each change of input }
    OffC : Array[1..6] of SpecArray; { Conc values for each change of input }

    Ind, Dep : integer;      { Main loop counters }
    i, j : integer;          { secondary counters }
    Flag : Boolean;          { error flag }
    HalfY : integer;
    o     : integer;         {offset counter}
    dummy : double;          { intermediate result }
    dum2  : double;
    DM2   : double;          {Dmult squared}

Const

    Offset : Array[1..6] of real = (-2.5,-1.5,-0.5,0.5,1.5,2.5);


begin
   {Initialize all vectors to Zero }
   For i := 1 to S.NComp do
     DeltaX[i] := 0;

   For j := 1 to S.NSpec do
     DeltaC[j] := 0;

   HalfY := NumY div 2;
   Predict := False;
   DM2 := Dmult * Dmult;

   InitArrays(TRUE);

  { Now Calculate the derivative for each change in XT or LogX, and sum}
  { the squares in DeltaX and DeltaC }

  Ind := 1;                       { Calculate the squared derivative for }
                                  { X and C for each TotX or LogX change }
    Repeat
      For o := 1 to 6 do begin
        InitArrays(FALSE);
        If (S.comps[ind].CType = CTotDiss) then
           Tot[Ind] := Tot[Ind] + S.Comps[Ind].StdDev * Dmult * Offset[o]
          Else LogX[Ind] := LogX[Ind] + S.Comps[Ind].StdDev * Dmult *
                            Offset[o] * LogConvr;

         Flag := SolidSolve;
         If (Flag = False) then Exit;

         OffX[o] := LogX;
         OffC[o] := Conc;
       end;                              {end of for statement}

      Case Order of
       1: begin
         For i := 1 to S.NComp do
          DeltaX[i] := DeltaX[i] + (OffX[4,i]-OffX[3,i])*(OffX[4,i]-OffX[3,i])
                       /(Dm2);
         For j := 1 to S.NSpec do
          DeltaC[j] := DeltaC[j] + (OffC[4,j]-OffC[3,j])*(OffC[4,j]-OffC[3,j])
                       /(Dm2);
         end;
        2: begin
          For i := 1 to S.NComp do begin
            dummy := OffX[5,i] + 6*OffX[4,i] - 6*OffX[3,i] - OffX[2,i];
            DeltaX[i] := DeltaX[i] + dummy * dummy / (81*DM2);
           end;
          For j := 1 to S.NSpec do begin
            dummy := OffC[5,j] + 6*OffC[4,j] - 6*OffC[3,j] - OffC[2,j];
            DeltaC[j] := DeltaC[j] + dummy * dummy / (81*Dm2);
           end;
          end;
      3: begin
           For i := 1 to S.NComp do begin
             dummy := OffX[6,i]+7*OffX[5,i]+12*OffX[4,i]-12*OffX[3,i]
                     -7*OffX[2,i]-OffX[1,i];
             dum2 := OffX[6,i] + OffX[1,i] - OffX[4,i] - OffX[3,i];
             DeltaX[i] := DeltaX[i] + dummy * dummy / (1444*DM2)
                          + dum2 * dum2 / (32 * DM2 * DM2);
           end;
           For j := 1 to S.NSpec do begin
             dummy := OffC[6,j]+7*OffC[5,j]+12*OffC[4,j]-12*OffC[3,j]
                     -7*OffC[2,j]-OffC[1,j];
             DeltaC[j] := DeltaC[j] + dummy * dummy / (1444*Dm2);
          end;
         end;
         end;                            {end of case statement }


     Ind := Ind + 1;

  Until(Ind > S.NComp);

  Ind := 1;                 { Calculate the squared derivative in X and C}
                            { For each change in Log K }
  Repeat
      For o := 1 to 4 do begin
        InitArrays(False);
        LogKZero[Ind] := LogKZero[Ind] + S.Specs[Ind].StdDev * Dmult
                         * Offset[o] * LogConvr;

        Flag := SolidSolve;
        If (Flag = False) then Exit;

        For i := 1 to S.NComp do
          OffX[o,i] := LogX[i];
        For j := 1 to S.NSpec do
          OffC[o,j] := Conc[j];

      end;                      { end of for }

      Case order of
      1: begin
          For i := 1 to S.NComp do
            DeltaX[i] := DeltaX[i] + (OffX[4,i]-OffX[3,i])*
                         (OffX[4,i]-OffX[3,i])/(Dm2);
          For j := 1 to S.NSpec do
            DeltaC[j] := DeltaC[j] + (OffC[4,j]-OffC[3,j])*
                         (OffC[4,j]-OffC[3,j])/(DM2);
         end;
      2: begin
          For i := 1 to S.NComp do begin
            dummy := OffX[5,i] + 6*OffX[4,i] - 6*OffX[3,i] - OffX[2,i];
            DeltaX[i] := DeltaX[i] + dummy * dummy / (49*DM2);
           end;
          For j := 1 to S.NSpec do begin
            dummy := OffC[5,j] + 6*OffC[4,j] - 6*OffC[3,j] - OffC[2,j];
            DeltaC[j] := DeltaC[j] + dummy * dummy / (81*DM2);
           end;
          end;
      3: begin
           For i := 1 to S.NComp do begin
             dummy := OffX[6,i]+7*OffX[5,i]+12*OffX[4,i]-12*OffX[3,i]
                     -7*OffX[2,i]-OffX[1,i];
             DeltaX[i] := DeltaX[i] + dummy * dummy / (1444*DM2);
           end;
           For j := 1 to S.NSpec do begin
             dummy := OffC[6,j]+7*OffC[5,j]+12*OffC[4,j]-12*OffC[3,j]
                     -7*OffC[2,j]-OffC[1,j];
             DeltaC[j] := DeltaC[j] + dummy * dummy / (1444*DM2);
          end;
         end;
        end;                                   {End of Case Statement}

     Ind := Ind + 1;

  Until(Ind > S.NSpec);

  For i := 1 to S.NComp do                     {Calculates SD's from}
    DeltaX[i] := Sqrt(DeltaX[i])/LogConvr;     {Sum square derivatives}
  For j := 1 to S.NSpec do
    DeltaC[j] := Sqrt(DeltaC[j]);

  {For i := 1 to HalfY do
   begin
    MCVars[i] := 1;
    If (Y[i] <= S.NSpec) then SDs[i] := DeltaC[Y[i]]
      else SDs[i] := DeltaX[Y[i]-S.NSpec];
   end;  }

  InitArrays(False);
  Predict := SolidSolve;
  {RecordResult(Means); }

{  ClrSCr;
  For i := 1 to HalfY do
     Writeln('Mean ',i,' = ',Means[i],'    Error = ',SDs[i]);
  Readln;}

end;


{*************************************************************************}
{      Combined Sweep with Monte Carlo simulation at each point           }
{      Exit Codes                                                         }
{         1- no data set present                                          }
{         2- no titration parameters                                      }
{         3                                                               }
{*************************************************************************}

Function MoSweep : integer;

var ttr : integer;                { incremental counter for titration points}
    i, com, spe   : integer;      { counters }
    ErrorCode : integer;          { error message }

Begin

  ErrorCode := 0;
  InitArrays(TRUE);                { Initialize ALL computational arrays }

  ttr := 0;                        { prepare to titrate }

  repeat                           { Main titration Loop }

    ttr := ttr + 1;                { increment counter  }

   With S do
    case TiType of                { change the X variable in the Data Set }

     TiTotal: Comps[TitrX].TotalX := XData[ttr];  { Class 1 component sweep }

     TiFree: Comps[TitrX].GuessX := XData[ttr];  { Class 2 component sweep }

    end;                          { end of cases for change X }
                                  { end of With S }


    If (MCreps = 0) then
     begin                      { If no Monte Carlo simulation desired }
       InitArrays(False);       { Now initialize all computational arrays, }
                                { retaining the previous FreeX concentrations }
                                { for all Known Total components }

   If (not SolidSolve) then Exit;  {if a calculation does not converge, Exit }
                                   {and return False}

    Replace;
    for com := 1 to S.NComp do     { record molarity of monitored comps, specs}
           YData[com,ttr] := S.Comps[com].Conc;
    for spe := 1 to S.NSpec do
           Ydata[spe + S.NComp, ttr] := S.Specs[spe].Conc;

   end
     {Else If (MCReps in [1..3]) then       {Error Analysis by Derivative}
     {  begin
         HalfY := NumY div 2;
         Flag := Predict(MCReps);
         For i := 1 to HalfY do
           begin
             Ydata[i,ttr] := Means[i];
             YData[i+HalfY, ttr] := SDs[i];
           end;
       end}
     Else begin                 { Monte Carlo Simulation }
           {HalfY := NumY div 2;

          { If SAVE then begin
              str(ttr, CountString);
              MonteName := MonteFileName + CountString;
            end;    }

           {ErrorCode := MonteCarlo;
           {FindMeans;
           For i := 1 to HalfY do
             begin
               YData[i, ttr] := Means[i];
               YData[i+HalfY, ttr] := SDs[i];
             end;                  }
        end;                       { End of Else }

  until (ttr = TitrPts);  { end titration loop }
  MoSweep := ErrorCode;
 end;                                        { End of Function MoTitrate }

{***********************************************************************}
{      Combined Titration with Monte Carlo simulation at each point       }
{      to save an error surface for one variable                          }
{       Exit Codes
{         1- no data set present                                          }
{         2- no titration parameters                                      }
{         3                                                               }
{*************************************************************************}

Function Surface : integer;

var ttr : integer;                { incremental counter for titration points}
    i   : integer;                { counter }
    flag  : integer;              { error message }
    HalfY : integer;
    BinValue : real;
    TF : Boolean;

Begin

 { if (Not SAVE) then
     begin
        Surface := 3;
        Exit;
      end;            }

  InitArrays(TRUE);                { Initialize ALL computational arrays }

  HiBin := 30;                     { Prepare to collect error distribution }
  LoBin := -30;

 { ClrScr;
  TF := TopGetReal('Lower Bound for Distribution',LoBin, Hibin, Lobin);
  TF := TopGetReal('High Bound for Distribution',Lobin, Hibin, Hibin);
  NumBins := TopGetInt('Number of distribution bins',10,MaxBins,NumBins);
  Range := (Hibin-Lobin)/NumBins;     }

  MonteName := MonteFileBase;
  Assign(MonteFile, MonteName);    { Open file to store results }
  Rewrite(MonteFile);

  Write(MonteFile,'Bins  X  Zmean ZSD');
  For i := 1 to NumBins do begin
    write(MonteFile,' Z',i);
   end;
  Writeln(MonteFile);

  ttr := 0;                        { prepare to titrate }

  repeat                           { Main titration Loop }

    ttr := ttr + 1;                { increment counter  }

   With S do
    case TiType of                { change the X variable in the Data Set }

     TiTotal: Comps[TitrX].TotalX := XData[ttr];  { Class 1 component sweep }

     TiFree: Comps[TitrX].GuessX := XData[ttr];  { Class 2 component sweep }

    end;                          { end of cases for change X }
                                  { end of With S }



                                   { Monte Carlo Simulation }
      HalfY := 1;
      {Flag := MonteSurface;
      FindMeans;
      YData[1, ttr] := Means[1];
      YData[2, ttr] := SDs[1];

      Write(MonteFile,Lobin+ttr*range:12,'  '); {Col1 = Bin values}
      {Write(MonteFile,XData[ttr]:12,'  ');      {Col2 = titration variable}
      {Write(MonteFile,'  ',Means[1],'  ',SDs[1]); {Col3 and Col4 = mean & SD}
      {For i := 1 to NumBins do
         Write(MonteFile,'  ',Distrib[i]);      { then save distribution }
      {Writeln(MonteFile);      }

  until ((ttr = TitrPts) or (Flag > 0));  { end titration loop }

  For i := (TitrPts+1) to NumBins do
     Writeln(MonteFile, Lobin+i*BinWidth:12);
  Close(MonteFile);

  Surface := Flag;
 end;                                        { End of Function Surface }



end.
